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Two basic correlation functions are calculated for a model of N harmonically interacting identical 
particles in a parabolic potential well. The density and the pair correlation function of the model are 
investigated for the boson case. The dependence of these static response properties on the complete 
range of the temperature and of the number of particles is obtained. The calculation technique 
is based on the path integral approach of symmetrized density matrices for identical particles in a 
parabolic confining well. 
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I. INTRODUCTION 

Generalizing the Feynman approach of identical particles in a box |ffl] to the case of identical particles in a parabolic 
confining potential, the present authors derived analytic expressions for the propagator and for the partition function 
of a system of N harmonically interacting identical particles (bosons or fermions) in a parabolic well 0], hereafter 
referred to as I. This model, giving rise to repetitive Gaussian integrals, also allows to obtain the generating function 
for the partition function. For an ideal gas of non-interacting particles in a parabolic well, this generating function 
coincides with the grand-canonical partition function. For interacting particles this generating function circumvents 
the constraints on the summation over the cycles of the permutation group at the expense of doing an extra path 
integral. 

In the present paper the one-point and two-point correlation functions of the model are calculated using their 
generating function as we did for the thermodynamic properties of the model. Also here we have to introduce extra 
path integrals of Gaussian nature to facilitate the cyclic summations. 

The one-point and two-point correlation functions of the model are obtained for the boson case as well as for the 
fermion case. But in view of the recent interest in Bose-Einstein condensation in a trap |3| |^], the explicit evaluation 
and the discussion of the results are restricted to the boson case in the present paper. The fermion case will be studied 
in a forthcoming paper. 

In the case of distinguishable particles, the correlation functions play a key role in the variational approximation for 
path integrals This variational method can be reformulated for indistinguishable particles, and the knowledge of 
the one- and two-point correlation functions for harmonic trial actions is as crucial as it is for distinguishable particles. 
For any algorithmic approach to many-body diffusion for interacting particles, the knowledge of the correlation 

functions of the model is very useful to test the actual implementations. Furthermore, the model provides an example 
of an exactly tractable system with interactions, which clearly exhibits the effects of Bose-Einstein condensation in 
the specific heat [[y] and in the moment of inertia 

The one-body potential energy V± and the two-body potential energy V2 of the model system are given by 



V = V 1 + V 2 ; * = ^X> 2 ; ^-^E^-r,) 2 . (1.1) 

i=i 3,1=1 



The two-body interaction is assumed to be repulsive; replacing — J 1 by oj 2 in V2 gives the attractive case. As a result 
of the diagonalization one obtains in each dimension one degree of freedom (the center of mass) with frequency f2, 
and N — 1 degrees of freedom with frequency w given by 



w 



y/Q 2 - Ncj 2 . (1.2) 
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For this many-body system distinguishability of the particles therefore implies that one is dealing with a system that 
reduces to 3iV degrees of freedom, each degree of freedom representing one linear harmonic oscillator. It is clear that 
for such a system the propagator, the thermodynamic functions and the correlation functions are well known 

For identical particles (bosons or fermions) the propagator can be obtained from the decomposition of the underlying 
processes in terms of 4 orthogonal processes with well-defined boundary conditions [pHUl. A typical sample path 
for fermions is provided by the subprocess with absorption at the boundary for the ^-direction (leading to a fermion 
diffusion process), while for the y- and the z-direction a boson diffusion process with reflection at the boundary has 
to be used. This procedure is used to generate the trajectory of the walkers (including the y- and z -components) . 
The path of a walker in this particular subprocess terminates if the motion in euclidean time gets absorbed along 
the rr-direction. Indistinguishability has therefore the important effect of making the coupled oscillator problem in 3 
dimensions also a genuine 3D problem. 

In I we showed how extra degrees of freedom for the center of mass can be introduced (in Fourier space) leading 
to a propagator which factorizes. The extra degrees of freedom can be integrated out afterwards. As a consequence, 
the internal degrees of freedom of the interacting oscillator system can be considered as the degrees of freedom of 
another non-interacting oscillator system. This mapping allows to use the grand-canonical partition function of the 
non-interacting system as the generating function for the system in interaction, provided the fugacity and hence the 
thermodynamical potential are identified as usual. 

The calculation of the density and the pair correlation function heavily relies on the calculations presented in I, 
which makes it difficult to make this paper self-contained without repeating some of the material presented in I. We 
tried to partly overcome this inconvenience by using the same notation as in I, and by a limited number of explicit 
references to the detailed manipulations in I if similar situations are encountered. 

The paper is organized as follows. In the next section the one- and two-point correlation functions are calculated 
for identical particles (bosons or fermions). In the third section the density and the pair correlation function are 
analyzed for the boson case. In the fourth section we discuss the results and draw some conclusions. 



II. STATIC RESPONSE PROPERTIES OF A MANY-BODY SYSTEM 

For the calculation of the static response properties of a many-body system, the correlation functions (e lq r! ) 7 
and Ylw ( e * q '( ri ~ ri ')) / are the key ingredients. The subscript I emphasizes that identical particles are considered 
(which can be specified to be bosons with subscript B or fermions with subscript F) in 3D. In the path integral 
approach the expectation values of an expression A (f ', r) are given by 

» J df J df'Kj (r, f}\r>, r) A (r' , r) K z (r>, r|r, 0) 
{A (r ' T))l = JdrKj (f,P\f,0) ' (2J) 

where Kj is the statistical propagator of the identical particles and r is the 37V-dimensional vector containing the 
coordinates i"i , • • • , r of the N particles. In this notation the probability density, the pair correlation function and 
their Fourier transforms are given by 

»(*) = ^(l>(*-*)) = / ^sV^-^ ^Ev^Oz' ( 2 - 2 ) 



* M = ^ ( E * ( r - r < + r„)\ = / T^W"^ — 9q = 1 £ . (2.3) 

Collecting the appropriate expressions for the propagators Kj (r, /3|f, r) and K] (r', r|r, 0) from I, one sees that the 
Fourier transforms n q and g q are given by 

N 



q NZ r J J (2ir) 3 



/dre-* ? ^e^-i ? ^C P n^(( J Pr) j , / 3|r,,0) , (2.4) 



l p j=i 



N N 
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where Zi is the partition function, K (r^-, /3 1 r j , 0) is the propagator of a 3D harmonic oscillator with frequency w, P 
denotes a permutation matrix and £ = — 1 assures the required anti-symmetry for fermions, whereas £ = +1 describes 
bosons. 

We first show how a tractable expression can be obtained for n q , and subsequently use an analogous procedure to 
calculate <? q . 



A. The single-particle expectation values 



Using the cyclic decomposition, and denoting by M/ the number of cycles of length £, the expression for n q can be 
written in terms of the cycles in the same way as we did for the partition function in I: 



1 



dRdk 

5"' 

(2^) 3 



ik-R 



Y Y £M e lC e (k, q) \,,„ M (AC/ (k))^" 1 TT -i jz- AC/, (k)) M *' 



(2.6) 



where 



/C,(k,q)= / dr t+1 drt- - f dv 1 5{T i+1 ~v 1 )e^X[K{T J+1 ^\v ] ,0) w e-^ r > 1 (2.7) 

3=1 

and AC/ (k) = AQ (k, q = 0) is precisely the same function as found in Eq. (2.20) of I for the determination of the 
partition function. Both the integrations over k and R only involve Gaussian integrands, and one eventually finds 
after some algebra 



n q = exp 



hq 2 ( coth|/37if2 coth h/3Hw 



4mTV V 



(2.8) 



The factor in front of n q accounts for the center of mass, and obviously reduces to unity in the non-interacting case 
where w — fl. The factor n q itself denotes the expectation value of J2i e lq r ' in the subspace of the relative coordinate 
system only with partition function Zj (TV) : 



1 



TVZ/ (TV) 



Mi — M A 



~y^ Mgl exp 



Hq 2 coth ^£(3hw 
Amw 



n 



{l-Y)M t 



1 



M e l£ Mi \2smhU/3hw 



3 M e 



(2.9) 



In AQ (k, q) one recognizes the partition function (over a time interval iff) of a driven 3D harmonic oscillator with 
frequency w 



AC, (k, q) = J dvK (r, ip\T, 0) w B- C *-«,(r).r(r) . fq (r) = ,i_ k £ tf (r _ + ^ ^ _ q W 



(2.10) 



which is known |l|,|6| in closed form: 

1 



3 q 27 7n 2mw 



f q (r) ■ f q (a) cosh (f - l T - CT l) ^ 



2sinh^) ' i! sinh±£/?to 



AC/ (k,q) 

The calculation of <£> q , given / q (r) as a sum of delta-functions, is straightforward. The result is 
and consequently: 



(2.11) 



h ( £k 2 1 k • q 1 1 

i-q — — -; ( coth-(3hw — 2—— coth -fihw + q 2 coth -£f3hw I , 

Amw \ N z 2 TV 2 2 



(2.12) 
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JC e (k, q) = JCe (k) exp 
JCe (k) ! 



/ ftk ■ q 1 hq 2 1 

— — coth -pnw — coth -£8hw 

\2Nmw T 4mw 2 H 



(2sinhf hwj 



■ exp 



ft ^ 2 , 1 _ , 



Introducing the generating function Q\ (u, q) = X)jv=o I^ 7 (-^0 ^^q] u 



N . 



5iKq)=E E 

jV=0 M X ---M N 



^] M e £ exp 



ftg 2 coth ^£/3Hw 
Amw 



11 /if J 



£(/-!), 



the summations can be done: 



£1 («, q) = Sj (it) 



Ll Me! L^(2sinhi£/3ftu;) c 
^ ^P (-^ coth \Z(3hw) 

—11 

(2sinh^/3ftw) 3 



M, 



(2.13) 
(2.14) 



(2.15) 



(2.16) 



where Sj (it) = X)iv=o ^/ (^0 * s ^ ne generating function of the partition function Z/ (AT) of TV identical oscillators 
in the relative coordinate system, studied in I. Consequently 



-y 



N 



exp 



hq 2 



coth MBKu 



h (N-£) 



N 



(2smh±lf3hw) 3 Z/W 
Considering the limit q — >0, it should be noted that the sum rule n q= o = 1 is indeed satisfied. 



(2.17) 



B. The two-particle expectation values 

Similarly as for the treatment of the single-particle correlation function, the Fourier transform which allows to treat 
the center-of-mass coordinate as an independent degree of freedom is introduced. The cyclic decomposition of the 
permutations implies that a factor e* qri occurs once in each position of each cycle. Furthermore a different factor 
gjq r,/ occurs m each position of each cycle which differs from r; (i.e. the case I = I' has to be excluded if r; and ri> 
are within the same cycle). Taking these bookkeeping considerations into account, the cyclic decomposition of the 
summation over the permutations leads to 



9q 



dlidk 



R e n 

Mi-Mjv \ e 



£(t-l)M e 

M e l£ Mt 



ZiN J J (2irY 

( E -ll Kt (k, q;j + 1) {K t (k))^" 1 Ue^e (k)) M *' 



E* M < 



\ 

+£ (Me - 1) K t (k, q) JCe (k, -q) (JC e (k)) Me ~ 2 Ue^e ( k )) Mf ' 
\+Y.v ¥ 4l'M t .K l (k,q)AC* (k,-q) (jC t (fc))^" 1 (JC e (k))^'" 1 1^,,^^, (£,« (k)) M "' ] 



(2.18) 



where JCe (k, q) and (k) are defined as in the previous subsection, and a function JCe (k, q;j) is introduced which is 
given by 



/Q(k,q;j) - / rfr £+1 / dr e ■ ■ ■ / dr 1( J (r /+1 - n) e^^e"^ J] ^ (r /+1 , /?|r y , 0) w e"^^' . 



(2.19) 



In (k, q;j) one recognizes the partition function (over a time interval £f3) of a driven 3D harmonic oscillator with 
frequency w 



e-i 



JCe (k, q;j + 1) = / drX (r, £d\r, 0) w e" ^ drh,W)-r(r) . hq (T; j} = f _L k £ 5 (r _ _ zq(5 (t) + iqS (r _ _ 

3 i'=o 

(2.20) 
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Similarly as for the single-particle correlation functions this expression is known MJ6J in closed form: 



Kt (k, q;j + 1) 



2 sinh ■# Tito 



e *,(i); * q (i) = |/ dr/ da 



h q (r,j)-h q (r,j) cosh b 



It — erj ] fi.u> 



2mw 



sinh M(3hw 



The explicit evaluation of the influence function VPq (j) is somewhat involved but straightforward: 

£ Ki^e^ + l hq 2 cosh \£f3hw ~ cosh (\l-j)l3hw 



*q CO = 



iV 2 4mro e^ fi ™ - 1 2mw 



sinh ^£(3Hw 



and consequently 



/C^ (k, qjj + 1) = Ki (k) exp - 



% 2 cosh — cosh — j) 

2mw sinh ^£(3hw 



(2.21) 



(2.22) 



(2.23) 



Using the results obtained in the previous subsection for ICe (k, q) and ICi (k) , the Fourier transform of the pair 
correlation function becomes: 



1 



■9q 



Z,N 



dRdk 

T' 

(2^) 3 



R x e n 

Mi—Ms \ t 



MfM M * 



JJ(/C, (k)) M * 



2^=1 ex P 



^g 2 cosh ^ifihw— cosh^^— jjfihw 
2mw sinh ±t(3hw 



+^ (M< - 1) exp (-| ^ coth yphw 

K + Ei'ju f ' M t' ex P (-jt£; ( coth + coth h e P hw )) ) 



(2.24) 



The condition ^^IMi — N on the cyclic decompositions simplifies \\i (k)) Mf , and the integrations over k and 



R are then straightforward, resulting in ( g^g^n^l 



. The remaining summation over the cycles can again be done 



if one introduces the appropriate generating function 

oo 

G 2 («,q)= $1 WN 9q \' 



(2.25) 



JV=0 



which lifts the restriction on the number of cycles of given length. The summation is straightforward. With b = e & w 
one obtains: 



G 2 (u, q) = 5/ (u) < ~Z /J exp 



Hq 2 (1 



3 = 1 



2mw 



Tig 2 1 + b 6 
Amw 1 — ¥ 



Using (XXi a if = Y1T=2 ELi °j a ^-i and defining 



Qtj (b) 



l-b e 



the terms can be combined into 



£2 (u, q) = S 7 (u) ^ — —3- Y 



exp 



(1 - 6») (1 - 6*-i) 
7m? 2 1 



2mw j (6) 



+ £(Q £j (&)) J exp(-^-Q £j (&) 



ft<7 2 



(2.26) 



(2.27) 



(2.28) 



It should be noticed that only cycles with length at least two contribute to the pair correlation function, as is to 
be expected. Because the series expansion of Q2 (u, q) in powers of u yields Zj (AT) N <? q as the coefficient of u N one 
obtains immediately 
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£ — 2 V / ^ — 1 

In the case of the pair correlation function g q the sum rule .g q =o = N — 1 can also be checked. The proof proceeds by 
induction, but it is rather tedious and is omitted here. All details on this calculation are provided upon request. 

III. BOSON DENSITY AND PAIR CORRELATION FUNCTION 

In this section the density and the pair correlation of the model are evaluated for the boson case, and the conden- 
sation effects on these quantities are studied. 

A. Density 

The density n (r) in the case of boson statistics can be obtained from n q and reads: 

1 AZ/(JV-£) fwm , \3/2 / mwr 2 \ . . 

with 

Ae ~ coth \t(3hw + i (ff coth pffi - coth pftw) ' ^'^ 
where r stands for the distance from the center of the confining potential. The density is ccntrosymmetric, as a 
consequence of the isotropy of the model. Introducing the activity pj as in I from Zb (N) = j-jw & ^ 1 _ b j^ IljLo Po tnc 
density can be rewritten as follows 



- > 1 / wmAi\ 3 ^ 2 ( raw , . \ (l — 6 J ) 3 

r>^— t(^T^) exp(-— r 2 ^) fl i ^- J "' ,:, - :,) 



1 \ - 1 / wmAi \ ( raw , , \ -i-r 



j=iV-M-l 

which allows for a recursively defined expression well suited for numerical evaluation 



,s 1(1-^) (i-b N - i y (i-b N - 2 Y h-b 2 ) 6 (i-bf ,,,, 

n r = T7 A oi + A k + A a 3 + + \a N -i + L a N 3.4 

N PJV \ PJV-l \ PN-2 V P2 \ Pi 



with 

where p = r^/ 2 ^ is a natural dimensionless quantity proportional to the distance from the center. Since T = 
tT c <~ tN 1 / 3 (where T c is the condensation temperature for the Bose-Einstein transition) p/N 1 / 6 is a natural quantity 
against which to plot the density. The results are summarized in two figures. In Fig. 1, the density n (0) /n,T=o (0) 
in the origin is shown (where tit = q (0) is the density in the origin at zero temperature), and exhibits a pronounced 
dependence on the condensation temperature. In Fig. 2, n (r) jn (0) is plotted as a function of r for 1000 particles. 
For comparison, the corresponding densities for the case of distinguishable particles are plotted in Fig. 3. T c is only 
used as a reference temperature for comparison purposes to Fig. 2; it does not have the meaning of a condensation 
temperature if the particles are distinguishable. 

The typical shape of the density as a function of the temperature is shown in Fig. 4 for T = 0, T = 0.9T C , T = T c 
and T = 1.1T C , where the spatial dependence of the density n(x,0, z) is plotted at a fixed value y = for 1000 
particles. It should be noted that the sudden appearance of an intense peak below T c when sweeping through the 
condensation temperature is manifestly present also in isotropic systems. 

The center-of-mass contribution to the density can be substantial for a limited number of particles. For 1000 
particles this single degree of freedom quantitatively has a negligible contribution to the density as a function of r^/w; 



the effects of the interaction enter in the eigenfrequency w = V^ 2 — Nlu 2 which determines the scaling parameters in 
the figures. 



G 



B. Pair correlation function 



An analogous analysis as for the density can be made for the pair correlation function: 

1 y^Zi(N-£) ^- 1 6# < ^ rmw „ \3/2 f / mwr 2 ^ \ _ ( mwr 2 1 \1 . , 

w U S^ h(-TT ««J + f (~5T ft")]' (3 ' 6) 

Introducing 

E 1 / mm \ 3 / 2 1 

(with q\ — 0) a recurrence relation for g (r) can be obtained similarly as for n (r) , but with at from ([T^) replaced by 
qi. In the same units as for the density, g (r) /<? (0) is plotted in Fig. 5. Similarly as for the density in the previous 
subsection, g (r) jg (0) of distinguishable particles is shown in Fig. 6 for comparison to the boson case in Fig. 5. 

The interpretation of Fig. 5 requires some caution, because g(r) / g (0) is plotted, and the magnitude of g (0) 
strongly depends on the condensation temperature. Nevertheless, it is clearly seen that the probability to find 
another particle at a relatively small distance r from some particle is very pronounced in the condensate. Above the 
critical temperature a more substantial contribution is obtained at relatively large distances, but the boson character 
still manifests itself by a larger probability to find particles at relatively short distances from the center. 

IV. DISCUSSION AND CONCLUSIONS 

The present paper concludes the boson part of a study of an exactly soluble model containing one type of particles. 
The system exhibits condensation at a finite temperature T c . The thermodynamical quantities such as the internal 
energy, the specific heat and the moment of inertia have been studied before by the present authors and also by 
others fl^-[lq| as far as some non-interacting aspects or ground state properties jl9| are concerned. The density of 
this model is an important response property. Precisely the concentration variations as a function of the cooling and 
the field are invoked to establish the condensation transition . Therefore it is comforting that the predictions of 
this theoretical model -which of course constitutes a simplification- bare some resemblance to the simulated density 
of an anisotropic boson oscillator model |po| and with the experimental situation in several aspects. It should be 
noted (i) that the magnetically induced anisotropy of the trap is not taken into account in the present paper, and (ii) 
that the interparticle interactions are replaced by harmonic two-body interactions. Furthermore it is assumed that 
the time scale of the experiment allows for an interpretation in terms of the thermal equilibrium response properties. 
It is clear that these simplifications deserve further investigations. 

Also the pair correlation function of the model is an important quantity especially if one wants to investigate the 
modifications due to a more realistic interparticle interaction using a variational approach. Indeed, an estimate of 
the effective interparticle potential along the lines of ((jjTj requires the pair correlation function of this exactly soluble 
model, that is taken as the zeroth order approximation to the system. 

From the methodological point of view, the projection and the generating function technique allows to obtain 
tractable and exact expressions for the thermodynamic quantities and for the static response functions of a system of 
identical boson oscillators. 
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Figure captions 

Fig. 1: Boson density n (0) /nr=o (0) a t the origin as a function of temperature. 

Fig. 2: Scaled density of bosons n (r) /n(0) for 1000 particles as a function of the distance r from the center for 
several temperatures. 

Fig. 3: Scaled density n (r) /n(0) for 1000 distinguishable particles for comparison to Fig. 2. 

Fig. 4: Probability density n (x, 0, z) for 1000 particles as a function of x and z for T = 0, T = 0.9T C , T — T c and 
T = 1.1T C . 

Fig. 5: Scaled pair correlation function g(r) /g (0) of 1000 bosons as a function of the distance r for several temper- 
atures. 

Fig. 6: Scaled pair correlation function g (r) / g (0) of 1000 distinguishable particles for comparison to Fig. 5. 
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